﻿namespace SpikingAnalyticsLib

open System.IO
open System.Collections.Generic
open SpikingNeuronLib
open GeneralUtilities
open IOUtilities
open PatternUtilities
open SpikingAnalyticsLib.PatternExtensions

module internal NetworkUtilities =

    // ===============================================================================================================
    // Network utility functions

    let RunNetwork (spn:NetworkBase) runSeconds (backgroundPattern:'b option) (inputPattern:IStimulus option) =
        let background = if backgroundPattern.IsSome then backgroundPattern.Value else null

        if inputPattern = None then
            spn.Run(runSeconds, background)
        else
            // with background, connection weights to non-stimulated neurons fall to zero (assuming STDPDepressionAmplitude > STDPPotentiationAmplitude)
            // without background, connection weights to non-stimulated neurons gradually increase by 0.01/sec as weight dervatives are always zero
            spn.Run(runSeconds, background, inputPattern.Value)

    // ===============================================================================================================
    // Connection utility functions (CrossbarNetwork)

    // e.g. to get outgoing connections i.e. connections keyed by pre-synaptic neuron
    // use: GetConnectionsByNeuron spn (fun connection -> connection.PreNeuron)
    // or: GetConnectionsByNeuron spn (fun connection -> connection.PostNeuron) for incoming connections
    let GetConnectionsByNeuron (spn:CrossbarNetwork) (selectNeuron:Connection -> int) =
        spn.Connections
        |> Seq.fold (fun (state:Map<int, Connection list>) connection ->
                let neuron = selectNeuron connection
                let connectionList =
                    if state.ContainsKey(neuron) then
                        connection :: state.[neuron]
                    else
                        [connection]
                Map.add neuron connectionList state
            ) Map.empty

    let GetPostNeurons neuron (outgoingConnections:Map<int, Connection list>) =
        outgoingConnections.[neuron]
        |> Seq.map (fun connection -> connection.PostNeuron)

    let ShowConnectedPostNeurons neurons (outgoingConnections:Map<int, Connection list>) =
        let PrintConnected neuron (outgoingConnections:Map<int, Connection list>) =
            printfn "Neurons connected to %d:" neuron
            let connectedNeurons =
                GetPostNeurons neuron outgoingConnections
                |> StringifyWithSeparator " "
            printfn "%s" connectedNeurons
        neurons
        |> Seq.iter (fun neuron -> PrintConnected neuron outgoingConnections)

    let GetConnectionsWithSharedPostNeurons preNeurons (incomingConnections:Map<int, Connection list>) (outgoingConnections:Map<int, Connection list>) =
        let GetConnectionsWithPreNeurons (preNeurons:Set<int>) postNeurons (incomingConnections:Map<int, Connection list>) =
            postNeurons
            |> Seq.map (fun neuron ->
                    incomingConnections.[neuron]
                    |> Seq.filter (fun connection -> preNeurons.Contains(connection.PreNeuron))
                )
            |> Seq.concat
                
        let sharedPostNeurons =
            preNeurons
            |> Seq.map (fun neuron ->
                    outgoingConnections.[neuron]
                    |> Seq.map (fun connection -> connection.PostNeuron)
                    |> Set.ofSeq
                )
            |> Set.intersectMany

        GetConnectionsWithPreNeurons (Set.ofSeq preNeurons) sharedPostNeurons incomingConnections

    let GetConnectionsDescription (connections:seq<Connection>) =
        connections
        |> Seq.map (fun connection -> connection.PreNeuron, connection.PostNeuron, connection.Delay)

    let PrintConnections (connections:seq<Connection>) =
        printfn "Pre\tPost\tDelay"
        GetConnectionsDescription connections
        |> Seq.iter (fun (preNeuron, postNeuron, delay) -> printfn "%d\t%d\t%d" preNeuron postNeuron delay)

    let GetPostNeuronsFromConnections (connections:seq<Connection>) =
        PrintConnections connections
        connections
        |> Seq.map (fun connection -> connection.PostNeuron)
        |> Seq.distinct
        |> Seq.toList

    let ShowConnections neurons (connectionMap:Map<int, Connection list>) =
        neurons
        |> Seq.iter (fun neuron ->
                connectionMap.[neuron]
                |> PrintConnections
            )

    // ===============================================================================================================
    // Connection utility functions (IzhikevichNetwork)

    // Create a map of post-synaptic neurons keyed by each pre-synaptic neuron
    let GetPostNeuronMap (network:IzhikevichNetwork) =
        let post = network.PostSynapticConnections
        let weights = network.SynapticWeights
        Seq.init (post.GetLength(0)) (fun preNeuron ->
                let postNeurons =
                    Seq.init (post.GetLength(1)) (fun columnIndex ->
                        let postNeuron = post.[preNeuron, columnIndex]
                        let weight = weights.[preNeuron, columnIndex]
                        let delay = network.GetDelay(preNeuron, postNeuron)

                        postNeuron, weight, delay
                    )
                    |> Seq.toList
                preNeuron, postNeurons
            )
        |> Map.ofSeq
        
    // ===============================================================================================================
    // Frame profiling utilities

    // Generate a map between selected neurons and their response histograms
    // where the response histogram is a map of spike counts keyed by temporal offset from an input pattern neuron
    let GenerateProfiles numTimeSlots patternStartNeuronID totalNeurons selectedNeurons firingEventList =
        let selectionMask = CreateFilter selectedNeurons totalNeurons

        firingEventList
        |> Seq.scan (fun state (time, nindex) ->                            // collect zero or more response spikes following each pattern spike
                if nindex = patternStartNeuronID then                       // new input
                    Some(time), None      // startTime, responseTime 
                else
                    let startTime, _ = state                                // new response
                    if startTime.IsSome then
                        startTime, Some(time - startTime.Value, nindex)     // return time difference (response - input) for each response neuron
                    else
                        None, None
            ) (None, None)
        |> Seq.choose (fun (_, response) -> response)                       // select responses
        |> Seq.filter (fun (time, nindex) -> selectionMask.[nindex])        // process only selected neurons
        |> Seq.fold (fun (dict:Map<int,int list>) (time, nindex) ->         // a list of reponse times for each neuron
                let times =
                    if dict.ContainsKey(nindex) then
                        dict.[nindex]
                    else
                        []
                Map.add nindex (time :: times) dict
            ) Map.empty
        |> Map.map (fun nindex times ->                                     // turn each list of response times into a separate map of counts keyed by time
                times
                |> Seq.fold (fun (dict:Map<int,int>) time ->
                        if dict.ContainsKey(time) then
                            let count = dict.[time]
                            Map.add time (count + 1) dict
                        else
                            printfn "Warning: response time (%d) was out of range (%d time slots allocated)" time dict.Count
                            dict
                    ) (PreAllocateMap numTimeSlots 0)
            )   

    // Produce a mapping: Map<neuron, (peakTime, peakCount)>
    // that quantifies the overall response to the input pattern
    let QuantifyResponseProfiles allMaps =
        let GetPeak responseTimes =
            let time, count =
                responseTimes
                |> Map.toSeq
                |> Seq.maxBy (fun (time, count) -> count)
            time, count
        allMaps
        |> Map.map (fun nindex responseTimes -> GetPeak responseTimes)

    // Generate all possible firing time combinations for 2 inputs (and therefore two firing times)
    let GetFiringTimeCombinations neurons maxDelay () =
        seq {
            for ft1 in 0..maxDelay-1 do
                for ft2 in 0..maxDelay-1 do
                    yield [ ft1; ft2; ]
        }
        |> Seq.map (fun firingTimeList ->
                Seq.zip firingTimeList neurons
            )

    // Generate all possible firing time combinations for 3 inputs (and therefore three firing times)
    let GetFiringTimeCombinations3 neurons maxDelay () =
        seq {
            for ft1 in 0..maxDelay-1 do
                for ft2 in 0..maxDelay-1 do
                    for ft3 in 0..maxDelay-1 do
                        yield [ ft1; ft2; ft3; ]
        }
        |> Seq.map (fun firingTimeList ->
                Seq.zip firingTimeList neurons
            )

    // Get the times from the specified sequence of firing events
    let GetEventTimes firingEvents =
        firingEvents
        |> Seq.map (fun (time, nindex) -> time)
        |> Seq.toList

    // Find the peak responses for all combinations of timing delays for the specified input neurons
    let GenerateResponseGrid (inputNeurons:int list) getCombinations spn runSeconds maxDelay inputPeriod selectedNeurons totalNeurons (collector:OneSecondTickDataCollector) =

        let patternGenerator = Some(new PoissonProcessPatternGenerator(totalNeurons, 1))
        let neuron1 = inputNeurons.[0]

        let responsePatterns =
            getCombinations ()
            |> Seq.mapi (fun index patternEvents ->
                    let eventTimes = GetEventTimes patternEvents
                    let eventTimesString = StringifyWithSeparator " " eventTimes
                    printfn "Testing input %d: %s" index eventTimesString

                    let inputPattern = Stimulus.Create(new Pattern(patternEvents), inputPeriod) :> IStimulus

                    collector.Clear()
                    RunNetwork spn runSeconds patternGenerator (Some(inputPattern))

                    let result = GenerateProfiles 250 neuron1 totalNeurons selectedNeurons collector.AllEventPairs

                    eventTimes, (QuantifyResponseProfiles result)
                )
            |> Map.ofSeq

        responsePatterns

    // Using a heuristic, determine which output neurons respond for all combinations of timing delays for a pair of input neurons
    // The heuristic uses low resolution (firing event) data to determine neuron responsiveness
    let GenerateLoResHeuristicResponseGrid neuron1 neuron2 spn runSeconds maxDelay inputPeriod selectedNeurons totalNeurons (collector:OneSecondTickDataCollector) =

        let GetRespondingNeurons frameSize minSpikeCount selectedNeurons totalNeurons firingEventList =
            let selectionMask = CreateFilter selectedNeurons totalNeurons
            let baseIndex =
                selectionMask
                |> Seq.findIndex (fun flag -> flag)

            let responses =
                firingEventList
                |> Seq.filter (fun (time, nindex, isBackground) -> selectionMask.[nindex])        // process only selected neurons
                |> Seq.map (fun (time, nindex, isBackground) -> time / frameSize, time % frameSize, nindex - baseIndex)       // remap time into each frame of size 'frameSize'; nindex -> 0-based index
                |> Seq.fold (fun (dict:Map<int * int,int>) (frame, time, nindex) ->         // key:nindex; value:number of spikes for each neuron
                        let spikes =
                            if dict.ContainsKey (frame, nindex) then
                                dict.[frame, nindex]
                            else
                                0
                        Map.add (frame, nindex) (spikes + 1) dict
                    ) Map.empty
                |> Map.filter (fun nindex spikeCount -> spikeCount >= minSpikeCount)

            responses

        let patternGenerator = Some(new PoissonProcessPatternGenerator(totalNeurons, 1))
        let getCombinations = GetFiringTimeCombinations [ neuron1; neuron2; ] maxDelay 

        let minSpikes = 2
        let frameSize = 1000    // one second frames
        let selectedFrame = 1   // skip frame zero

        let responsePatterns =
            getCombinations ()
            |> Seq.mapi (fun index patternEvents ->
                    let eventTimes = GetEventTimes patternEvents
                    let eventTimesString = StringifyWithSeparator " " eventTimes
                    printfn "Testing input %d: %s" index eventTimesString

                    let inputPattern = Stimulus.Create(new Pattern(patternEvents), inputPeriod) :> IStimulus

                    collector.Clear()
                    RunNetwork spn runSeconds patternGenerator (Some(inputPattern))

                    let responseMap = GetRespondingNeurons frameSize minSpikes selectedNeurons totalNeurons collector.AllEventTriplets
                    let respondingNeurons =
                        responseMap
                        |> Map.toSeq
                        |> Seq.filter (fun ((frame, nindex), _) -> frame = selectedFrame)
                        |> Seq.map (fun ((frame, nindex), spikeCount) -> nindex)
                        |> Seq.toList

                    if respondingNeurons.Length > 0 then
                        printfn " %d" respondingNeurons.Length
                    else
                        printfn ""

                    eventTimes, respondingNeurons
                )
            |> Map.ofSeq

        responsePatterns

    // Using a heuristic, determine which output neurons respond for all combinations of timing delays for a pair of input neurons
    // The heuristic uses high resolution data (membrane potential data (also weight data)) to determine neuron responsiveness
    let GenerateHiResHeuristicResponseGrid neuron1 neuron2 spn runSeconds maxDelay inputPeriod outputNeurons totalNeurons (collector:OneMillisecTickDataCollector) =

        // Return a map of (weightIndex, weightSaturated) pairs keyed by output neuron index
        let GetWeightSaturationKeyedByOutputNeuron threshold sampleCount layerOutputs (collector:OneMillisecTickDataCollector) =
            let TestIfWeightsSaturated threshold count (collector:OneMillisecTickDataCollector) =
                let startIndex = collector.GetNumberOfSamples("MembraneData").Value - count
                collector.ConnectionDataPairsRowView.Value
                |> Seq.map (fun (weightIndex, row) ->
                        let saturated =
                            row
                            |> Seq.skip startIndex                                      // last 'count' values in each row...
                            |> Seq.forall (fun (_, weight) -> weight >= threshold)    // ...all over threshold...
                        weightIndex, saturated
                    )
            // assumes ALL connections in this layer are selected
            let connectionsPerOutputNeuron =
                let layerConnections = collector.NumberOfSampledConnections
                layerConnections / layerOutputs
            TestIfWeightsSaturated threshold sampleCount collector
            |> Seq.groupBy (fun (index, saturated) -> index / connectionsPerOutputNeuron)      // key by output neuron index
            |> Seq.map (fun (nindex, responses) -> nindex, Seq.toList responses)
            |> Map.ofSeq


        // For each output neuron, test if the neuron spiked and with what spiking pattern
        // return a sequence of (neuronIndex, patternSpiked, spikePattern) triplets
        // where patternSpiked is true if the neuron spiked at any time during the sampling period
        // where spikePattern is the spike pattern during the sampling period but with reduced resolution
        let TestIfNeuronsSpiked threshold (collector:OneMillisecTickDataCollector) =
            let ReduceResolution factor pattern =
                let result =
                    pattern
                    |> Seq.mapi (fun index spiked -> index, spiked)             // add an index
                    |> Seq.groupBy (fun (index, spiked) -> index / factor)      // reduce resolution by factor
                    |> Seq.map (fun (index, window) ->                          // return true if any spikes in window
                            window
                            |> Seq.exists (fun (_, spiked) -> spiked)
                        )
                    |> Seq.toList
                result
            collector.MembraneDataPairsRowView.Value
            |> Seq.map (fun (neuronIndex, row) ->
                    let spikePattern =
                        row
                        |> Seq.map (fun (v, _) -> v >= threshold)
                        |> ReduceResolution 20
                    let patternSpiked =
                        spikePattern
                        |> Seq.exists (fun spiked -> spiked)
                    neuronIndex, patternSpiked, spikePattern
                )
            |> Seq.toList

        // Select the neuron indices that responded to the current input pattern
        // i.e. select output neurons where all weights have saturated
        // and there are at least 'minSpikes' spikes
        let SelectActiveNeurons (weightSaturationByNeuron:Map<int, _>) minSpikes data =
            data
            |> Seq.filter (fun (neuronIndex, _, pattern) ->
                    let allSaturated =
                        weightSaturationByNeuron.[neuronIndex]
                        |> Seq.forall (fun (weightIndex, saturated) -> saturated)
                    let MinimumSpikes pattern threshold =
                        let numSpikes =
                            pattern
                            |> Seq.filter (fun spiked -> spiked)
                            |> Seq.length
                        //printfn "%d" numSpikes
                        numSpikes >= threshold
                    allSaturated && (MinimumSpikes pattern minSpikes)
                )
            |> Seq.map (fun (neuronIndex, _, _) -> neuronIndex)
            |> Seq.toList

        let saturationThreshold = 9.5
        let sampleCount = 1000   // test weight saturation on the last 'sampleCount' samples
        let spikeThreshold = 30.0
        let minSpikes = 5

        let patternGenerator = Some(new PoissonProcessPatternGenerator(totalNeurons , 1))
        let getCombinations = GetFiringTimeCombinations [ neuron1; neuron2; ] maxDelay 

        let responsePatterns =
            getCombinations ()
            |> Seq.mapi (fun index patternEvents ->
                    let eventTimes = GetEventTimes patternEvents
                    let eventTimesString = StringifyWithSeparator " " eventTimes
                    printf "Testing input %d: %s" index eventTimesString

                    let inputPattern = Stimulus.Create(new Pattern(patternEvents), inputPeriod) :> IStimulus

                    collector.Clear()
                    RunNetwork spn runSeconds patternGenerator (Some(inputPattern))

                    // Note that neuron and connection indices are in terms of the data collector (i.e. the selected neurons and connectors rebased to zero)
                    let weightSaturationByNeuron = GetWeightSaturationKeyedByOutputNeuron saturationThreshold sampleCount outputNeurons collector
                    let respondingNeurons =
                        TestIfNeuronsSpiked spikeThreshold collector
                        |> SelectActiveNeurons weightSaturationByNeuron minSpikes

                    if respondingNeurons.Length > 0 then
                        printfn " %d" respondingNeurons.Length
                    else
                        printfn ""

                    eventTimes, respondingNeurons
                )
            |> Map.ofSeq

        responsePatterns

    // ===============================================================================================================
    // Data collection utility functions

    // Run the network with randomly-selected stimuli from inputPatterns
    // Return the stimuli selected for each frame
    let CollectFixedFrameDataRandomSelection (spn:CrossbarNetwork) (firingEvents:List<int * int * 'a>) runSeconds (backgroundPattern:'b option) (inputPatterns:Map<'c, IStimulus option>) =

        let background = if backgroundPattern.IsSome then backgroundPattern.Value else null
        let r = System.Random()
        let numPatterns = inputPatterns.Count

        // the inputPatterns map is allowed to have arbitrary keys
        // create a list of allowed pattern keys
        let patternKeys =
            inputPatterns
            |> Map.toSeq
            |> Seq.map (fun (key, value) -> key)
            |> Seq.toList

        let timeShiftedResults = new List<int * int * 'a>()

        let patternSelectionMap =
            Seq.init runSeconds (fun sec ->
                    let patternKey =
                        let selectedIndex = r.Next(numPatterns)
                        patternKeys.[selectedIndex]
                    let inputPattern =
                        inputPatterns.[patternKey]

                    firingEvents.Clear()
                    if inputPattern.IsSome then
                        spn.Run(1, background, inputPattern.Value)
                    else
                        spn.Run(1, background)

                    // copy time-shifted firing events
                    let timeShiftedFiringEvents =
                        firingEvents
                        |> Seq.map (fun (time, nindex, label) -> sec * 1000 + time, nindex, label)
                    timeShiftedResults.AddRange(timeShiftedFiringEvents)

                    printf "."  // show progress

                    sec, patternKey
                )
            |> Map.ofSeq

        patternSelectionMap, timeShiftedResults

    // Run the network with the specified stimuli and specified order
    let CollectFixedFrameDataOrderedSelection (spn:CrossbarNetwork) (firingEvents:List<int * int * 'a>) (backgroundPattern:'b option) (inputPatterns:Map<'c, IStimulus option>) (stimulusOrder:'c list) secondsPerPattern =

        let background = if backgroundPattern.IsSome then backgroundPattern.Value else null
        let numPatterns = inputPatterns.Count

        let patternSelectionMap =
            stimulusOrder
            |> Seq.mapi (fun index patternName ->
                    Seq.init secondsPerPattern (fun sec -> (index * secondsPerPattern) + sec, patternName)
                )
            |> Seq.concat
            |> Map.ofSeq

        let timeShiftedResults = new List<int * int * 'a>()

        let runSeconds =
            stimulusOrder
            |> Seq.fold (fun runSeconds patternName ->

                    let inputPattern = inputPatterns.[patternName]

                    firingEvents.Clear()
                    if inputPattern.IsSome then
                        spn.Run(secondsPerPattern, background, inputPattern.Value)
                    else
                        spn.Run(secondsPerPattern, background)

                    // copy time-shifted firing events
                    let timeShiftedFiringEvents =
                        firingEvents
                        |> Seq.map (fun (time, nindex, label) -> runSeconds * 1000 + time, nindex, label)
                    timeShiftedResults.AddRange(timeShiftedFiringEvents)

                    printf "."  // show progress

                    runSeconds + secondsPerPattern
                ) 0

        patternSelectionMap, timeShiftedResults

    // Run the network with randomly-selected patterns from inputPatterns
    // The input pattern is applied at a randomly chosen time in each one second run in the range 0 to 800 msec
    // Return the input patterns selected for each one second run
    let CollectFrameFreeRunData (spn:CrossbarNetwork) (firingEvents:List<int * int * 'a>) runSeconds (backgroundPattern:'b option) (inputPatterns:Map<int, Pattern option>) =
        // shift the specified pattern forward or backwards in time
        let ShiftPattern (inputPattern:Pattern option) shiftMilliSeconds =
            if inputPattern.IsSome then
                let shiftedEventArray =
                    inputPattern.Value.EventList
                    |> Seq.map (fun evt -> FiringEvent(evt.Time + shiftMilliSeconds, evt.NeuronIndex, evt.Label))
                    |> Seq.toArray
                Some(new Stimulus(1, shiftedEventArray))
            else
                None

        let background = if backgroundPattern.IsSome then backgroundPattern.Value else null
        let r = System.Random()
        let numPatterns = inputPatterns.Count
        let timeShiftedResults = new List<int * int * 'a>()

        let patternSelectionMap =
            Seq.init runSeconds (fun sec ->
                    let selectedIndex = r.Next(numPatterns)
                    let shiftMilliSeconds = r.Next(800)
                    let inputPattern = ShiftPattern inputPatterns.[selectedIndex] shiftMilliSeconds

                    firingEvents.Clear()
                    if inputPattern.IsSome then
                        spn.Run(1, background, inputPattern.Value)
                    else
                        spn.Run(1, background)

                    // copy time-shifted firing events
                    let timeShiftedFiringEvents =
                        firingEvents
                        |> Seq.map (fun (time, nindex, label) -> sec * 1000 + time, nindex, label)
                    timeShiftedResults.AddRange(timeShiftedFiringEvents)

                    sec, (selectedIndex, shiftMilliSeconds)
                )
            |> Map.ofSeq

        patternSelectionMap, timeShiftedResults

    // ===============================================================================================================
    // Network structure utility functions

    // Given a group of related firing events, produce a list of linked firing events
    // where each link is of the form [ (preTime, preNeuron); (postTime, postNeuron); connectionWeight ]
    let GetGroupLinksFromPostNeuronMap jitter (postMap:Map<int, (int * float * int) list>) groupFiringEvents =

        // create a list of firing times and associated probabilities for each neuron
        let neuronFiringTimeProbabilityMap =
            groupFiringEvents
            |> Seq.fold (fun (state:Map<int, (int * float) list>) (time, nindex, probability) ->
                    let timeList =
                        if state.ContainsKey(nindex) then
                            state.[nindex]
                        else
                            []
                    Map.add nindex ((time, probability) :: timeList) state
                ) Map.empty

        let linkedEventsList =
            groupFiringEvents
            |> Seq.map (fun (time, nindex, _) ->
                    postMap.[nindex]
                    |> Seq.filter (fun (postNeuron, postWeight, postDelay) -> neuronFiringTimeProbabilityMap.ContainsKey postNeuron)
                    |> Seq.map (fun (postNeuron, postWeight, postDelay) ->
                            // create a link to each firing event of the post neuron
                            // with the restriction that post firing events must come after the current firing event
                            // and the time difference must be restricted to some reasonable delay: axonal delay + integration time
                            let postFiringTimes =
                                neuronFiringTimeProbabilityMap.[postNeuron]
                                |> Seq.filter (fun (postFiringTime, _) ->
                                    let isPostEvent = postFiringTime > time
                                    let isWithinJitter =
                                        let firingTimeDifference = postFiringTime - time
                                        (firingTimeDifference >= postDelay) && (firingTimeDifference <= postDelay + jitter)
                                    isPostEvent && isWithinJitter
                                    //(isPostEvent) && ((postFiringTime - time) < maxDelay)
                                )
                            postFiringTimes
                            |> Seq.map (fun (postFiringTime, postProbability) -> ((time, nindex), (postFiringTime, postNeuron), postWeight, postProbability))
                        )
                    |> Seq.concat
                )
            |> Seq.concat
            |> Seq.toList

        let minTime, maxTime =
            let timeList =
                groupFiringEvents
                |> Seq.map (fun (time, _, _) -> time)
                |> Seq.toList
            Seq.min timeList, Seq.max timeList

        minTime, maxTime, linkedEventsList

    // Using a list of outgoing connections for each neuron and a list of firing events
    // deduce a connectivity graph
    let GetGroupLinksFromPostConnectionMap jitter synapticThreshold (outgoingConnections:Map<int, Connection list>) groupFiringEvents =

        // get a list of connected post neurons for each pre neuron plus the connection weight between the two neurons
        // connection strengths must be above the specified synaptic threshold (typically they need only be non-zero)
        let postNeuronMap =
            outgoingConnections
            |> Map.map (fun preNeuron connectionList ->
                    connectionList
                    |> Seq.filter (fun connection -> connection.Weight > synapticThreshold)
                    |> Seq.map (fun connection -> (connection.PostNeuron, connection.Weight, int connection.Delay))
                    |> Seq.toList
                )

        GetGroupLinksFromPostNeuronMap jitter postNeuronMap groupFiringEvents

    let ParseLinkedEvents(lines:seq<string>) =
        let data =
            lines
            |> Seq.map (fun line -> Tokenize (List.ofArray (line.ToCharArray())))
            |> Seq.toList

        let firstFiringTime, lastFiringTime = ConstructTuple<int * int>(data.[0])
        let linkedEvents =
            data
            |> Seq.skip 1
            |> Seq.map (fun (tokenList:Token list) ->
                    if tokenList.Length = 4 then
                        let t1, n1, t2, n2 = ConstructTuple<int * int * int * int>(tokenList)
                        (t1, n1), (t2, n2), 0.0
                    else if tokenList.Length = 5 then
                        let t1, n1, t2, n2, w = ConstructTuple<int * int * int * int * float>(tokenList)
                        (t1, n1), (t2, n2), w
                    else
                        failwith "Expecting a row of four integers, or four integers and a float"
                )
            //|> Seq.cast<SpikingAnalyticsLib.LinkedEventPair>
            |> Seq.map (fun ((t1, n1), (t2, n2), w) ->
                let preEvent = new FiringEvent(t1, n1, EventLabel.Foreground)
                let postEvent = new FiringEvent(t2, n2, EventLabel.Foreground)
                new LinkedEventPair(preEvent, postEvent, w, 1.0)
            )
            |> Seq.toList

        firstFiringTime, lastFiringTime, linkedEvents

    // Analyse the network connectivity to find a set of neurons with shared output neurons
    // For each set of neurons found, use the axonal delays to determine an appropriate input pattern
    let FindCongruentPatterns numPatterns (network:CrossbarNetwork) =

        let FindCongruentPattern (random:System.Random) outputNeuron maxDelay network =

            // work backwards: start with an output neuron and find shared inputs
            let MakeLayers (random:System.Random) outputNeuron (network:CrossbarNetwork) =

                let rec MakeLayers (random:System.Random) layer maxLayers numExcitatoryNeurons (incomingConnections:Map<int, Connection list>) outputLayerNeurons (excludeSet:Set<int>) (delayMap:Map<int, int>) =

                    let inputConnectionMap =
                        let SelectConnections numInputs (incomingConnectionMap:Map<int, Connection list>) =
                            // select numInputs values for each output neuron key
                            incomingConnectionMap
                            |> Map.map (fun output connections ->
                                    let rec SelectIndices numConnections numRequired (selectedIndices:Set<_>) =
                                        if selectedIndices.Count = numRequired then
                                            selectedIndices
                                        else
                                            let selected = random.Next(numConnections)
                                            SelectIndices numConnections numRequired (selectedIndices.Add selected)

                                    let incomingExcitatoryConnections =
                                        connections
                                        |> Seq.filter (fun connection -> connection.PreNeuron < numExcitatoryNeurons)
                                        |> Seq.toList

                                    // randomly select the required number of connections
                                    let len = incomingExcitatoryConnections.Length
                                    if len >= numInputs then
                                        let selectedIndices = SelectIndices len numInputs Set.empty
                                        selectedIndices
                                        |> Seq.map (fun index -> incomingExcitatoryConnections.[index])
                                        |> Seq.toList
                                    else
                                        // Insufficient connections (pattern will be discarded)
                                        incomingExcitatoryConnections
                                )

                        outputLayerNeurons
                        |> Seq.map (fun neuron ->
                                let incomingConnectionList =
                                    incomingConnections.[neuron]
                                    |> Seq.filter (fun connection -> not (excludeSet.Contains connection.PreNeuron))
                                    |> Seq.toList
                                neuron, incomingConnectionList
                            )
                        |> Map.ofSeq
                        |> SelectConnections 3

                    printfn "Layer %d" layer
                    inputConnectionMap
                    |> Map.iter (fun outputNeuron incomingConnectionList ->
                            let connectionDelays =
                                incomingConnectionList
                                |> Seq.map (fun connection -> connection.PreNeuron, connection.Delay)
                                |> StringifyWithSeparator " "
                            printfn "%d: %s" outputNeuron connectionDelays
                        )

                    if layer = maxLayers then
                        inputConnectionMap, delayMap
                    else
                        let preNeuronSet =
                            inputConnectionMap
                            |> Map.toSeq
                            |> Seq.map (fun (output, connections) -> connections)
                            |> Seq.concat
                            |> Seq.map (fun connection -> connection.PreNeuron)
                            |> Set.ofSeq

                        let delaysFromNeuronMap =
                                inputConnectionMap
                                |> Map.toSeq
                                |> Seq.map (fun (output, connections) ->
                                        let downstreamDelay = delayMap.[output]
                                        connections
                                        |> Seq.map (fun connection -> connection.PreNeuron, int(connection.Delay) + downstreamDelay)
                                    )
                                |> Seq.concat
                                |> Map.ofSeq

                        MakeLayers random (layer + 1) maxLayers numExcitatoryNeurons incomingConnections preNeuronSet (Set.union excludeSet preNeuronSet) delaysFromNeuronMap

                let numExcitatoryNeurons = network.ExcitatoryNeurons
                let incomingConnections = GetConnectionsByNeuron network (fun connection -> connection.PostNeuron)
                let delaysFromNeuronMap = Map.add outputNeuron 0 Map.empty
                MakeLayers random 1 2 numExcitatoryNeurons incomingConnections (Set.singleton outputNeuron) (Set.singleton outputNeuron) delaysFromNeuronMap

            let inputConnectionMap, delaysFromNeuronMap = MakeLayers random outputNeuron network

            let inputPattern =
                let eventDelays =
                    inputConnectionMap
                    |> Map.toSeq
                    |> Seq.map (fun (output, connections) ->
                            let downstreamDelay = delaysFromNeuronMap.[output]
                            connections
                            |> Seq.map (fun connection ->
                                    let time = maxDelay - int(connection.Delay)
                                    let neuron = connection.PreNeuron
                                    downstreamDelay, time, neuron
                                )
                        )
                    |> Seq.concat
                    |> Seq.toList
                let maxDelay, _, _ =
                    eventDelays
                    |> Seq.maxBy (fun (downstreamDelay, time, neuron) -> downstreamDelay)
                eventDelays
                |> Seq.map (fun (downstreamDelay, time, neuron) -> (maxDelay - downstreamDelay) + time, neuron)
                |> Seq.toList

            inputPattern

        let rec RunTilSuccess (random:System.Random) outputNeuron totalNeurons maxDelay network =
            try
                let pattern = FindCongruentPattern random outputNeuron maxDelay network
                if pattern.Length = 9 then
                    printfn ""
                    Some(pattern)
                else
                    printfn "Pattern discarded: insufficient incoming excitatory connections\n"
                    let outputNeuron = random.Next(totalNeurons)
                    RunTilSuccess random outputNeuron totalNeurons maxDelay network                    
            with
                | :? System.InvalidOperationException ->
                    // unable to select 3 incoming connections (after filtering)
                    // try again with a different output neuron
                    let outputNeuron = random.Next(totalNeurons)
                    RunTilSuccess random outputNeuron totalNeurons maxDelay network

        let r = System.Random()
        let totalNeurons = network.TotalNeurons
        let maxDelay = network.MaximumAxonalDelay
        Seq.init numPatterns (fun index -> r.Next(totalNeurons))                // choose random output neuron
        |> Seq.map (fun outputNeuron -> RunTilSuccess r outputNeuron totalNeurons maxDelay network)
        |> Seq.toList

    // ===============================================================================================================
    // Miscellaneous utility functions

    let CreateKeyCountMap (data:seq<'Key * 'Value>) =
        data
        |> Seq.fold (fun (state:Map<'Key, int>) (key, value) ->
                let count =
                    if state.ContainsKey(key) then
                        state.[key]
                    else
                        0
                Map.add key (count + 1) state
            ) Map.empty

    // ===============================================================================================================
